District flood vulnerability assessment using analytic hierarchy process (AHP) with historical flood events in Bhutan

Flood hazards are common in Bhutan as a result of torrential rainfall. Historical flooding events also point to flooding during the main monsoon season of the year, which has had a huge impact in many parts of the country. To account for climate change patterns in flood hazards in Bhutan, 116 historical flood events between 1968 and 2020 for 20 districts were retrieved and reviewed. The preliminary review revealed that the frequency of flood occurrence has increased by three times in recent years. In this study, seven flood vulnerability (FV) indicators were considered. Five are the attributes of historical floods, classified into a number of incidents for flood events, fatalities, affected population, and infrastructure damages including economic losses. Additionally, the highest annual rainfall and existence of a flood map were other two indicators considered. Using historical data, flood hazard and impact zonation were performed. The analytic hierarchy process (AHP) method was employed to derive a multi-criteria decision model. This resulted in priority ranking of the seven FV indicators, broadly classified into social, physical/economic, and environmental. Thereafter, an indicator-based weighted method was used to develop the district flood vulnerability index (DFVI) map of Bhutan. The DFVI map should help researchers understand the flood vulnerability scenarios in Bhutan and use these to mediate flood hazard and risk management. According to the study, FVI is very high in Chhukha, Punakha, Sarpang, and Trashigang districts, and the index ranges between 0.75 to 1.0.


Introduction
Climate change's extreme events have triggered various natural hazards in Bhutan, causing huge impacts every year. Floods among other natural hazards continue to impend due to the effect of climate change [1]. Globally, flood hazards are one of the most destructive and costliest natural disasters [2][3][4][5]. Although flood damage is widespread and losses are enormous, due attention has yet to be paid in terms of hazard, vulnerability, and risk in Bhutan. The need for such a study has become important due to the limited number of studies conducted. Harm to property, deterioration of the environment, damage to infrastructure, and loss of life are some a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 of the consequences of flood hazards [6]. Also, the impacts of floods are increasing as the result of cascading climate change scenarios, with predicted global flood exposure to escalate by three times by 2050 [7]. The combination of torrential rainfall and topography, land use and land-cover pattern, and demographic status further substantiates Bhutan to be highly vulnerable to flood hazards.
The potential impacts of flooding to exposed elements such as human beings, communities, properties, and the environment are assessed through flood vulnerability. According to the World Bank [8], a vulnerability assessment is conducted to understand how systems such as infrastructures could be damaged or destroyed, how service and business-chain supply could be interrupted and how the community could suffer through property loss, negative health impacts, and fatalities. Vulnerability has been well defined by Devkota [9] and according to this, a system's predisposition to suffer damage during an extreme event is a representation of vulnerability and is therefore considered to be the magnitude of the impact on exposed system.
Based on local observations and drawing lessons from other countries, for Bhutan it is important to collate and analyse all the information linked to adverse effects of extreme weather phenomena over the past few decades that have led to many natural disasters, especially floods. The conceptual framework of flood vulnerability can be stated under various settings based on the purpose and perception of researchers [7]. In this study, the historical event-based framework was formulated to ascertain the hazard and impact levels in order to deduce and quantify vulnerability at the district level in terms of district flood vulnerability index (DFVI). The studies conducted by Diakakis [10] employed historical event-based hydrological data for flood history analysis and its contribution to hazard, while Burger et al. [11] accounted for climatic data. Discharge of past flood events was used by Sudhaus et al. [12], and the flow model for different return periods by Baky et al. [13], however, in some cases, historical events were often used only to enhance the knowledge or reconstruct the hydrological extremes [14][15][16]. In particular, historical traits are rarely mapped to assess the vulnerability of a region due to the absence of proper records. To this end, the current framework accounts for vulnerability driven by attributes of historical flood events that correlate with social, physical, economic, and environmental factors. The current study considered these factors as key indicators to derive flood vulnerability. An indicator-based approach has been widely used by many researchers and has been successfully employed in many regions at different scales, e.g., flood vulnerability indices at varying spatial scales [17], reducing the complexity of the flood vulnerability index [18], using the flood vulnerability index as a knowledge base for flood risk in urban areas [19], using the vulnerability index for urban flooding [20] and creating a new flood vulnerability index adapted for the pre-Saharan region [21]. Some also used a combination of flood vulnerability indicators and climate change impacts [22,23]. Additionally, the literature review indicates a similar framework used at the provincial level in China using multidimensional flood vulnerability. The framework is based on data envelopment analysis, which considered indicators such as population, fatalities, agriculture, and economic vulnerability [24]. The current study presents the historical flood-event records between 1968 and 2020 that include flash floods of both river floods and surface water floods, and glacial lake outburst floods (GLOF). The historical sequence of flood events not only provides a comprehensive listing, but also highlights impacts and signals the prevalence from likely probable flood hazard scenario to heavy rain disasters [25]. By evaluating historical flood records, the frequency of flood-hazard events and subsequent impacts were mapped and vulnerability assessments were conducted to develop flood vulnerability indices at the district level using AHP and an open source QGIS software. According to Pathak et al. [7], the quantification of vulnerability is meant to build a more resilient community and help policy-makers facilitate prospective directions to frame out interventions effectively, serve as a pivotal tool for risk management and identify suitable climate change adaptation strategies. The research also attempts to address several UN sustainable development goals (SDGs), e.g., good health and well-being, building resilient infrastructure, having sustainable cities and communities, and combatting climate change [26].
Most of the districts in Bhutan are regularly hit by floods that cause tremendous loss of property and infrastructures, including loss of life. The mapping of flood vulnerability is important for handling flood problems. The main objective of this paper is to prepare a DFVI map for Bhutan. The analytic hierarchy process (AHP) has been employed to derive a multicriteria decision model to determine the significance and influence of flood vulnerability indicators. The DFVI will provide insight into flood vulnerability and further attention can be paid to conduct in-depth study for the districts which are highly vulnerable to flood hazards. Additionally, a DFVI map will assist planners in adopting strategies at the district level for resource allocation, mitigation measures and enhancing resilience to flood hazards.

Study area
Bhutan is situated in the Himalayan region, which extends between 26˚45'N and 28˚10'N latitudes and 88˚45'E and 92˚10'E longitudes, with three major zones of the southern foothills, inner Himalayas and higher Himalayas. The altitude ranges from 1000 m to more than 6000 m from sea level, as shown in Fig 1. According to the Population and Housing Census 2017, Bhutan's population was projected to be 735,553 under 20 administrative districts and 205 local administrative units with 163,001 households [27]. The study area extends up to 38,394 km 2 with a population density of 19.5 as of 2020 [28]. The study area is characterised by four major river basins, the Punatshangchu, Wangchu, Manaschu, and Amochu basins. According to the historical events, most of the regions at sub-catchment zones are affected by floods except for the Punatshangchu basin, which was hit by GLOF three times affecting downstream corridors of the basin such as Punakha, Wangdiphodrang, and Sunkosh (Tsirang). The impact was recorded highest in Punakha. The vulnerability exists mainly due to inhabitants that populate many pocketed rural areas in most of the sub-catchment zones and have a relatively low coping capacity compared to the urban area. Further, this population's vulnerability is grossly upscaled with the additional infrastructure and service requirements, which are sheer attributes that come from livelihood needs. Fig 2 represents the demographic status of Bhutan that accounts for social and economic indicators. In the context of the proposed research, the demographic scenarios show substantial relative exposure to flood hazard risk (and other natural hazards). For instance, the rural population accounts for 62.2% of the total population, while other building categories such as masonry, mud and timber structures, and temporary houses, are as high as 67.7%. Also, 34.7% of the population has not attended schooling, and 36.7% is economically inactive apart from significant contributions through other attributes. These socio-economic attributes highlight Bhutan to be highly vulnerable to flood hazards, and susceptible to impending extremes of climate change.

Historical flood data
The flood data were collected primarily through the National Centre for Hydrology and Meteorology (NCHM), the district administration office, and the archive of the national newspaper, Kuensel. Some of the data were also retrieved from respective web portals. From the data collection process, 116 historical annual flood events that occurred between 1968 and 2020 were retrieved ( Table 1). The data retrieved were de-clustered as per the administrative boundaries (districts) and classified into the number of incidents for flood events, fatalities, affected population, and infrastructure damages. The economic loss was classified based on the extent of damages. Since the historical records do not have complete details of loss assessment, the cumulative extent of the damage level for the respective district was considered to classify the loss (very low to very high).

Rainfall data
The precipitation in Bhutan manifests in monsoon characteristics and reveals significant variation, with the lowest total rainfall in January and the highest in July for most parts of the country [29]. According to Bhutan's climatic record from 1996, the precipitation trend has been increasingly pronounced over the years, resulting in numerous extreme events [30][31][32][33]. From  Bhutan State of the Climate 2020 [33], the highest rainfall records were used for all the districts, which served as an input parameter in QGIS to develop the rainfall map. The inverse distance weighting (IDW) method was employed for interpolating at random locations of the study area to show the spatial variations in five classes (Fig 4). The annual accumulated rainfall of 2020 is observed to be higher than the cumulative average annual rainfall between 1996 and 2019.
The southern belt of the country, which includes Sarpang, Samtse, Chhukha, and Samdrupjongkhar, usually received much higher rainfall, with moderate rainfall towards eastern and northern parts. The western and the central regions such as Bumthang, Paro, Punakha, Trongsa, and Wangdiphodrang experienced comparatively less rainfall. The annual cumulative rainfall is as high as 7220.30 mm in Sarpang and as low as 472.70 mm in Paro. In the extreme north of the greater Himalayan, Gasa has been also receiving moderately high rainfall annually. In 2020, Gasa recorded annual rainfall of 2508.30 mm.

Development of the flood vulnerability index
To articulate FVI, an indicator-based weighted method was employed, which is widely used for flood system vulnerability assessment in many regions across the globe, e.g., a assessment of urban vulnerability towards floods using the indicator-based approach in Santiago de Chile [34], flood impact in the Mekong Delta, Vietnam [35], an indexing approach to assess flood vulnerability in coastal cities of Mazandaran, Iran [36], assessment of coastal sites in the UK for flood vulnerability with a combined physical and economic index [37], social vulnerability assessment for flood risk analysis [38], and a physical vulnerability assessment on flash floods using an indicator-based methodology [39] and followed by application of multi-criteria decision analysis (MCDA), with the AHP model adopted by many researchers to integrate the priority of indicators, e.g., [40][41][42][43]. DFVI in the present study was conceptually derived from Nasiri et al. [44]. Alternatively, geospatial analysis and remote-sensing (RS) techniques are widely used for flood vulnerability assessments, considering different indicators under respective regional settings. Such studies are possible due to the availability of geographical information system (GIS) platforms and RS products, e.g., a remote sensing and GIS-based flood vulnerability assessment in Jiangxi province, China [45], urban flood vulnerability using AHP and GIS [46], social vulnerability assessment of flood risk using GIS-based MCDA [47,48], assessment of urban vulnerability to pluvial flooding using GIS applications [48], flood risk mapping using GIS and multi-criteria analysis [49], a GIS-based flood vulnerability assessment in Pasir Mas, Kelantan [50], a vulnerability assessment for flash floods using GIS spatial modelling in El-Arish city, North Sinai, Egypt [51], a flood vulnerability assessment using a GISbased multi-criteria approach in Attica region [52], and a coastal flooding risk assessment using a GIS-based spatial MCDA approach [53]. The outline of the method employed in the current study can be divided into four parts as follows: The historical data provided us with the attributes of flood hazards and impacts, which were broadly categorised under the three domains of social, physical/economic, and environmental factors, which served as the main criteria. Corresponding sub-criteria were identified and classified into five classes. These account for exposure and susceptibility to flood vulnerability. The existence of a flood map specifies the resilience indicator for a particular district. The FV indicators are as shown in Table 2.
After de-clustering the historical flood data records, the score was assigned to the corresponding class on an ascending order of scale from 1-5. Such a scale of the severity level of hazard or impacts was also used by Blistanova et al. [54], and many others elsewhere. Initially, this scale system was defined by Mchaughlin and Cooper [55], with '5' contributing most strongly to vulnerability, and '1' contributing least. The results of classification are subsequently utilised in developing hazard and impact maps at the district level of Bhutan. FV indicators were used in AHP to define pair-wise comparison to obtain priority indices which were then applied to computed DFVI. The study consists of seven FV indicators used in the AHP model under the umbrella of three base criteria. Details of criteria, sub-criteria, weighting, and the methodological process are presented in Fig 5. A multi-objective AHP developed by Saaty [56] employs multi-criteria decision analysis that uses a pair-wise comparison process to efficiently evaluate complex decisions by constructing the judgement matrix with the scale of absolute numbers 1-9 (Table 3). AHP uses hierarchical structures to represent a problem and then develops priorities for alternatives based on the judgement of the user [57]. The AHP model involves constructing a pair-wise matrix, eigenvalue and weighting coefficient calculation, and allows a check for consistency of the priority ranking by calculation of consistency ratio (CR) [58]. The consistency indices (CI) and CR of a given choice are calculated using Eqs 1 and 2.
where λ max is the maximum eigenvalue of the pair-wise comparison vector and n is the number of attributes. The random index (RI) is as indicated in Table 4.

Flood hazard zonation
The record of major events accounts for 116 historical floods in Bhutan between 1968 and 2020. The district-wide distribution of historical flood hazards is shown in Fig 6. The map also depicts the locations of the major flood events including recurrent events. Since 1990, the Chhukha district has been hit by the highest number of 16 flood events. A major flood disaster occurred in 1996 at Pasakha, when the Barsachu river flooded. The flood damaged more than 25 residential buildings, incurring losses of more than Nu. 30 million. Similarly, Phuentsholing city under the same district was also devasted by the 2000 Dhuti Khola flood and the 2016 Amochu flood. The frequency of floods in southern regions of the country shows a similar trend. Sarpang and Samtse experienced 10 and 15 major flood events respectively. The eastern part of the country also recorded major flood events, with 9 in Lhuentse, 10 in Trashiyangtse, and 12 in Trashigang. Other districts observed moderately few flood events, except for the Punakha district, which was hit by three major GLOF events. According to a study conducted by Gurung et al. [59], Punakha remains highly susceptible to GLOF due to active Lemthang Tsho and its association with climate change.

Flood hazard impacts
An assessment of four flood impacts (fatalities, infrastructure damages, affected population, and loss) was undertaken based on historical flood events (Fig 7). Such an impact assessment not only gives a detailed insight into flood characteristics but is also helpful in providing information on warning and evacuation [60,61]. To assess the flood hazard impact levels, the study considered cumulative impacts for all the FV indicators. Fig 7a shows the number of fatalities within each district due to flooding. The Chhukha, Punakha, and Sarpang districts recorded a greater number of fatalities, followed by five districts in the east. It should be noted that the greater number of fatalities is a consequence of a higher number of flood events, including major floods due to GLOF in Punakha. The central and extreme northern regions recorded a comparatively lesser number of flood events and fatalities. The Paro and Thimphu districts in the west recorded moderately higher fatalities. Infrastructure damage includes the total number of damage incidences of residential, commercial, industrial, and monumental buildings, retaining structures, water supply systems and tanks, instrumentation stations, bridges, roads, and irrigation channel segments. As shown in Fig 7b, the infrastructure damages were highest in the Sarpang and Trashigang districts, followed by a significant number in 10 districts. Eight other districts recorded a minimal number of infrastructure damages between 1 and 10.
Bhutan's population is densely populated, especially in urban areas (main district) under various sub-catchments. The total affected population was assessed based on the location of flood events in each district and the results are shown in Fig 7c. The assessment results indicate a lesser number of populations being affected in fewer districts. This shows that the affected populations are dependent not only on the frequency of flood occurrence but on the collateral impacts. For instance, the damage or collapse of a bridge could result in road network disruption, and the population in nearby local administrative units would thus be affected.
In the absence of detailed information, qualitative judgement was employed to classify economic losses (Fig 7d). The records of such information broadly describes damages and does not have complete assessment records. For example, loss of household properties or belongings, damages to agricultural farmland (wetland), dry land, cash crops, cattle, fruit trees, industrial setup, etc. These losses are often difficult to consider for an accurate valuation in terms of monetary value and that the damages have occurred in different time frames. Usually, a field assessment is a must to accurately estimate the actual loss. Hence, we rated the loss based on the overall damage extent from very low to very high in five classes. However, since an indicator-based weighted method was applied, the class and score approximately accede to these economic losses. The actual damage loss model accounts for the combination of direct tangible and indirect tangible costs including intangible costs, where effects are felt by the society, but the accompanying losses and damages are difficult to value in terms of monetary cost [62], something that is also highlighted and discussed in [63,64]. Here, economic losses are considered similar to intangible costs.

Flood vulnerability index
The indicator-based weighted method employs expert judgement to form pair-wise comparisons in AHP model execution. In the current study, pairwise priority was employed based on both expert judgement and the score of the class, derived from records of historical events. Table 5 represents the FV indicator scores for the Chhukha district that are used in the AHP pair-comparison ranking matrix, with the highest weighting of 5 in four criteria and a weighting of 4 for infrastructure damage. This also indicates that the FV indicators in the Chhukha district are one of the highest contributing attributes of historical flood events. These scores are correlated to pair-wise comparison for the AHP model to assign significance of importance with a scale of 1, 3, 5, 7, and 9 respectively, while a scale of 2, 4, 6, and 8 indicates values of intermediate importance. Conversely, inverse comparisons of less important variables were assigned 1 to 1/9. The priority matrix for the Chhukha district is shown in Table 6 to illustrate the importance of each FV indicator to overall flood vulnerability. The pair comparison matrix was formulated with 21 comparisons and a similar procedure was applied in 19 other districts. As a result of the AHP method, the district-level priority ranking for seven FV indicators is presented in Table 7. Each district displayed a unique correlation among the FV indicators. The economic loss, affected population and annual rainfall dominates the flood vulnerability scenarios, followed by infrastructure damage, and fatalities on account of the frequency of flood events. Some flood events indicate destructive nature of floods against frequency of occurrence. For example, Haa, Paro, Punakha, Thimphu, Wangdiphodrang, and Zhemgang are some of the districts which experienced such flood events. On the contrary, some of the other districts were impacted due to a higher number of flood events and rainfall intensity, e.g., Chhukha, Sarpang, Samtse Trashigang, and Trashiyangtse.
The priority ranking matrix requires a consistency ratio (CR) of less than 10%. In the few instances where the CR was more than 10%, pair-wise comparison indicators having the same score were re-prioritised based on the significance of the flood event and expert judgement. This re-prioritisation was observed to have a close correlation among the same group at the individual district level. The CR in AHP model analysis for 20 districts (Table 7) was achieved between 1.60% and 9.8%.
The representative normalised flood indicators for the Chhukha district are presented in Table 8. For existence of a flood map, 1 (Yes) was indicated, and 0 (No) was indicated for districts which do not have a flood map. The existence of a flood map is considered as a resilience indicator; however, the priority differed in each district based on significance of other indicators.  ------1   ------2   ------3   ---50 to 100  --4 30 to 50 20, The consolidated priority ranking of the AHP model for 20 districts in Fig 8 shows the unique correlation among indicators under each district, for example, in Samtse and Samdrupjongkhar districts, the environmental indicators of rainfall and frequency of flood account for~58%, imposing high susceptibility, compared to social and physical exposure indicators, which hold~16 to 25% respectively. On the contrary, in Sarpang and Chhukha, with a~25 to 37% contribution by all indicators, both the susceptibility and exposure levels have relatively similar impacts and these districts remain highly vulnerable to flood hazards. Also, Punakha displays a high social vulnerability of 46.9% and the overall vulnerability remains at an all-time high, as it falls in the regions of the largest Punatshangchu river basin with the potential for GLOF outbreaks. The Samtse district also recorded a greater number of floods with environmental indicators showing up to 58.4%, with flood vulnerability as high. However, the social exposure and impact are comparatively lower. Bumthang, Pemagatshel, Paro, and Wangdiphodrang remain high to flood vulnerability in terms of social and physical exposure, however, the overall vulnerability remains at low due to lesser rainfall pattern and flood events. In such a case, the consequence of a major flood event could be devastating. The rest of the districts remain at an all-time low to flood hazard with few seasonal impacts. The DFVI of each district was computed based on the general flood vulnerability index (FVI) formula (Eq 3), where, E, S, and R indicate exposure, susceptibility, and resilience respectively. This base formula is further re-generalised in Eq 4 to account for the number of indicators involved under each criterion, and DFVI is obtained by combining all the indicators using Eq 5 for each district. The results are presented in Table 9. Finally, the aggregated DFVI for each district is obtained by taking the sum of social, physical/economic, and environmental criteria using Eq 6 and normalising it by the maximum aggregated sum [11,44] as supplied in Table 10.
where W i is the weighted score and X i represents priority indices for corresponding indicators  respectively; The final obtained DFVI is scaled to define the vulnerability level [21] as presented in Table 11. The scaled vulnerability index provides flood vulnerability designation for each of the districts ranging between 0-0.01 as very low vulnerability to flooding, and 0.75-1.0 as very high vulnerability to flooding, with other intermediate values. The study indicates four districts with very high flood vulnerability, including two in the southern region, and one each in the east and west. One of the districts falls in the category of high vulnerability, another four in the medium category, and 11 districts in the low-flood vulnerability category. Districts with higher DFVI in the east, west and southern region of Bhutan remain constantly vulnerable to floods compared to other districts that are primarily contributed by susceptibility to a high frequency of events and heavy rainfall. As a result of the study, the DFVI map of Bhutan is presented in Fig 9.

Conclusions
In this paper, an assessment of historical flood events was undertaken and seven flood vulnerability indicators were identified under three broad categories of social, physical/economic, and environmental attributes. The multi-criteria decision analysis was carried out using AHP to formulate the influence of FV indicators to flood vulnerability in each district. The flood vulnerability index was computed based on an indicator-based weighted method to obtain the district flood vulnerability map of Bhutan. The following concluding remarks and recommendations are presented.
• According to the study, attributes of flood history show a unique function of susceptibility and exposure to flood vulnerability that play an important role in predicting the flood vulnerability index.
• Social indicators dominated the flood vulnerability in Chhukha, Punakha, Sarpang, and Trashigang, with 0.75-1.0 FVI, followed by Lhuentse with 0.57 FVI due to the high number of fatalities and the affected population. The Paro, Thimphu, and Trashiyangtse districts are also moderately high to social vulnerability to flood, with FVIs ranging between 0.25-0.50, while in the Samtse district, environmental indicators indicate a higher margin of flood vulnerability due to heavy rainfall and the significant number of flood events. The Bumthang, Dagana, Gasa, Haa, Mongar, Pemagatshel, Samdrupjongkhar, Trongsa, Tsirang, Wangdiphodrang, and Zhemgang districts remain at an all-time low to flood vulnerability. These districts receive comparatively low rainfall, with few flood histories. However, the study shows that exposure to flood vulnerability cannot be ruled out in these regions and FVI ranges between 0-0.1. The districts with a high vulnerability index indicate a very high possibility of reoccurrence of such flood events associated with extreme climatic conditions that are cascading every year.
• Social indicators significantly account for 15.8% to 46.9% of flood vulnerability. An in-depth study on social vulnerability of districts with high FVI is highly recommended to accurately assess the exposure and susceptibility level to flood hazards in Bhutan.
• The FVI developed in this research study can be used as a tool to strategise district-wise flood risk management and preparedness plans by the relevant stakeholders.